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Abstract 

We numerically study the evolution of a boost-invariant N = 4 SYM medium using AdS/CFT. 
We consider a toy model for the collision of gravitational shock waves, finding that the energy 
density first increases, reaches a maximum and then starts to decrease, matching hydrodynamics 
for late times. For the initial conditions we consider, the hydrodynamic scale governing the late 
time behaviour is to very good approximation determined by the area of the black hole horizon at 
initial times. Our results provide a toy model for the early time evolution of the bulk system in 
heavy- ion collisions at RHIC and the LHC. 
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I. INTRODUCTION 



The problem of colliding gravitational shock waves in spaces that are asymptotically 
Anti-de-Sitter has been of recent interest because it can serve as a toy model of the collision 
of two nuclei approaching at very high speeds. Hence it may provide — via the AdS/CFT 
conjecture l[ EJ — some qualitative insight in phenomena found in heavy-ion experiments 
at the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC). 

Boosting a mass at rest to very high velocities, its energy-momentum tensor in coordinates 
x ± = ^= becomes that of a gravitational shock wave, e.g. T ++ oc fx5(x + ), where /x is the 
energy per unit area which for a nucleus of atomic number a, radius R and Lorentz boost 



factor 7 = is 

am p ay/s NN 
/i0C7 ^ = ^^' (L1) 
where m p is the proton mass. It is well known how to model such an energy-momentum 

tensor using the AdS/CFT correspondence, namely by a metric of the form 3] (cf. jj-6|) 

2 -2dx + dx~ + (j)(z)5(x + )dx +2 + dx\ + dz 2 

US — - (l.ZI 
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where <f>(z) = — and we have set the AdS radius to unity. Here the description assumes 
the conjectured duality between N — 4 SYM at large coupling and large number of colours 
N c and classical gravity on AdS 5 x S 5 . Since Af = 4 SYM is a gauge theory, it behaves 
qualitatively similar to QCD, so some aspects of this work may translate to qualitatively 
similar phenomena found in nature. The constant k is usually set to but we will treat 
k as a free parameter to be adjusted at will in order to obtain a model that more closely 
resembles QCD. 

A collision of two nuclei can be modelled by a superposition of two shock waves in AdS§, 
moving in x + and x~ direction, respectively. While the line element before the collision is a 
simple superposition of the individual shock waves, the metric after the collision is in general 
hard to find. Unfortunately, while exact analytic solutions in four dimensional asymptotic 
Minkowski time have been derived many years ago |7| , no such solutions are known for shock 



waves of the form (11. 2D. There 



pioneered in 0, S], see also 
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ore, one has to resort to numerical techniques, which were 



111- 



Another aspect of Eq. (11.21) is that the collision of two such shock waves can be shown 
to violate Bjorken's conjectured invariance under rapidity boosts (|s|). Since experimental 



data for heavy-ion collisions does not seem to back up this invariance either, this can be 
considered a feature rather than a shortcoming of the present model, but at the price that 
the gravitational dynamics is 2+1 dimensional (rapidity, AdS radius and time) 10| rather 
than 1+1 dimensional. 

However, as we shall point out in the present work, it turns out that in the limit of weak 
shock waves ^<1, the leading order dynamics is in fact boost-invariant (cf. Ref. 12| ). The 
late-time behavior of such a strongly coupled boost-invariant M = 4 SYM medium has been 
known up to the 3 rd order in large r expansion 13Nl6|. However, one has to use numerical 



methods to fully understand the early-time properties of the system [171 ]. In this paper 
we use algorithms similar to those in Refs. p, 9] to solve Einstein's equations numerically 
in this approximation, and follow the evolution of the boundary energy-momentum tensor 
from the far-from equilibrium situation at early times to the hydrodynamic behavior at late 
times. Unlike Ref. [8|, |9j, we do not deform the boundary four dimensional metric of the 
AdS§ space but connect initial conditions derived analytically from the shock waves before 
the collision to the late time hydrodynamic regime. Our findings validate those of Ref. 111 ], 
where the authors use a different algorithm and start with arbitrary initial conditions. 

This paper is organized as follows. In Sec. Ullwe construct an ansatz metric function based 
on the approximate metric functions in the collisions of two weak shock waves. In Sec. II III 
two algorithms for numerically solving Einstein's equations in the bulk of the AdS§ space 
are described in detail. Our numerical results and an application to RHIC and LHC are 
presented in Sees. llVfVl In the Appendix [A] provide near-boundary power series expansions 
needed in our numerical calculations. 



II. COLLISIONS OF TWO WEAK SHOCK WAVES 

The line element (11.21) is highly singular at x = 0, and it is useful to first change 
to so-called Rosen coordinates x + = u , x~ — v + \4>(z)6(u) + | (<p'(z)) 2 u6(u) 2 , z = 
z + ^(p\z)u6{u) , with the result 

^ 2 -2dudv+dx 2 ± + (l+±cf>" (z)u8(u)Y dz 2 _ ^ 

S ~ (z+^'(z)ue(u)) 2 ' ' ' ' 

The collision of two shocks can be set up by superposing the above line element for one 
shock with an equivalent one for the other shock. The difficult part of the calculation then 



involves finding the line element in the forward light-cone. Using the standard matching 
conditions (metric needs to be continuous and piece-wise differentiable it has been pos- 



sible to find the metric in the approximation of small strength /x js]. Using the coordinates 
proper time f = \J1uv and space-time rapidity fj — | In - the result is given by {^J 

2 _ - [1 + K(f, fj, z)] df 2 + [1 + L(f, fj, z)] f 2 df] 2 + [1 + H(f, fj, z)] dx± 2 

z 2 [1 + 2z 2 fxf cosh(Y — fj)] 

| [1 + M(f, fj, z)] [1 + 6z 2 fif cosh(Y - fj)] 2 dz 2 
z 2 [1 + 2z 2 /ifcosh(y — fj)] 2 

where p, = /i/k and K,L,H,M were determined to be 

K(f, fj, z) = c x /i W - A ~z 2 + 0(f) 

L(f, fj, z) = Z^tf^ W - 5 -±^^~z 2 + 0{f) 

H(f, fj, z) = -2tfr 2 t - ^±^W + 0(f) 

o 

M(f, fj, z) = 16/Z 2 f 2 i 4 + W + Q{f) (2.3) 



by solving Einstein's equations R pv — jg^ u R ~ Gg^u = 0. Here c\ is a freely choosable 
integration constant that corresponds to some unfixed diffeomorphism freedom. Following 
Refs. , our numerical setup requires a line-element in the Eddington-Finkelstein form 
and, therefore, we have to transform to new coordinates r, rj, z. The relation between the 
old and new coordinates in the small ft limit may be found to be 

f = r + z - cosh \Y-rj\- ^iu 2 z\r + zf Cl + f cosh [2(Y - ??)]z 6 (r + z) 
^ z\175t 3 + 469r 2 z + 378rz 2 + 52z 3 ) - ^ Z g! + G (f) 



210 v ' ' ' 8(r + z 

/iz 4 sinh [Y - d /i 2 (8r 2 + 16z + lz 2 ) , , r rtr <ri , , 

77 = 77 — — — — L — — - + P v - V sinh [2(Y - 77)] + 0(/i 3 



5 = 2- 2/x2 3 (r + 2) cosh [y — 77] — ^^ 3(T + Z)4C1 + /x 2 cosh [2(Y - r])]6z 5 (r + z) 2 

6 

u 2 

-!—z 3 (25r 4 + 100r 3 z + 25rV - UQrz 3 - 119z 4 ) + 0(fi 3 ) . (2.4) 
30 

It turns out that to this order in the shock strength //, the metric in Eddington-Finkelstein 



coordinates is independent of 77 and hence boost-invariant in the sense of Bjorken 18| . Higher 



order corrections turn out to spoil this invariance, but it seems that — at least for weak 



shocks with /i < 1 — the initial dynamics is predominantly boost-invariant. In Eddington- 
Finkelstein coordinates, the line element can be parametrized in the following form 



ds 2 = 2drdr - Adr 2 + Y?e B dx\ + £ V 2B dr/ 2 , (2.5) 
where using r = - the metric functions are given as 

5r 4 3r 3 3r 2 + U ^ > ' 

B = -'log (h±IL\ + (612 + 7rr(234 + rr(217 + 75rr)))p' + 3) 

3 \ r J 315r 6 (l + rr) 

S 3 = r2(1+rT)+ P2 + lW5 rT )y +0( . 3) (26) 

A. Einstein's equations in Eddington-Finkelstein coordinates 

In Eddington-Finkelstein coordinates (12. 5p . Einstein's equations become 8] 

= ££' + 2£'£-2£ 2 , (2.7a) 

= EB'+|(E'S + B'E), (2.7b) 

= A // + 35 / J B - 12S'S/S 2 + 4, (2.7c) 

= S + i(5 2 S-A / S), (2.7d) 

= S" + |5 /2 £ , (2.7e) 

where for any function h(r, r) we have defined 

h = d T h + lAd r h, (2.8) 

and h' = d r h. Under the coordinate transformation 

r ->■ f — r — f, (2.9) 

one has 

A = A(r+ f, r) - 2^, £ = E(f + /, r), and B = B(r + f, r), (2.10) 

where / is an arbitrary function of r. It is easy to check that Einstein's equations are 
form-invariant under the above diffeomorphism. 



In the following, we will take ( 12.7dj) and (12. 7e[) as constraint equations and numerically 
solve fl2.7ap - fl2.7cp . which can be rewritten in the following form 

& = S, (2.11a) 

0' = 3B'S~h, (2.11b) 

S = Q6, (2.11c) 

B = -S~^<j), (2.11d) 

A" = 86S'S- 2 + 3B'(f)S-^ - 4, (2.11e) 

where S = S 3 . 



B. Apparent Horizons and Area of Trapped Surface 

The area of the trapped surface formed in the collision of two shock waves is of consid- 
erable interest since at late times, when the system is close to equilibrium, it can be used 
to extract the entropy of the system. Far from equilibrium, its physical interpretation is 
difficult [l6| but it is nevertheless interesting to track the area spanned by the apparent 
horizon, which is the location where out-going null vectors vanish. It should be pointed out 
that the apparent horizon is a local concept, coordinate-time dependent, and not invariant 
under coordinate transformations (space-time slicings). 

The location of apparent horizon may be calculated as follows: first determine the in 
and out-going null vectors l~, l + (corresponding to "light-rays" in ordinary space-time) from 
the condition g^l^lv = 0. Then find the apparent horizon from the criterion of vanishing 
expansion of the out-going null vectors, 

h ab V a lt = 0, (2.12) 

where h ab is the projected metric that is given by h ab = g ab — la lb . (The projected metric 
fulfills the requirement that multiplying l + by an arbitrary function B does not change the 
result f l2~T2|) V 



1. Before the collision 



Just before the collision of the two shock waves, where the line element is given by a 
superpositions of line elements of the form (12. ip . the location of the apparent horizon may 

6 



be calculated along the lines of [19]: parametrizing the surface at u = by v = —ipi(z), the 
normals to this surface are given by du, dv + dip\, so normal vectors Z M can be parametrized 
as l^dx^ = c\du + c 2 {dv + dipi) . The condition g^l^ly = at u = leads to the conditions 
c i = c 2 or C2 = , where the prime here denotes a derivative with respect to z. As a 

consequence we obtain a set of (out-going and in-going) null vectors normal to the surface 
at u = 

#=(^,1,0^®) , z; = (i,o,o ± ,o). 

(The constant function multiplying these vectors is arbitrary and has been set to unity.) 
Vanishing expansion implies 



□ 







which has the solution ipi(z) = \<t>(z) + c, where the constant c is unimportant for the 
following. To obtain the location of the trapped surface z = zh, we use the following 
boundary condition: we could have equally well started with the other shock wave and 
a surface at v — parametrized as u = —ip 2 {z)- Since at u = v = both surface normal 
vectors have to coincide, one finds ipi = ip 2 and iP'±(zh) = 2 = j(P' 2 (zh) and as a consequence 

~z H =(2tfY 1/ \ 

Now the "area" of the trapped surface is given by 

Ah = J \/ det 9ab\ s dx ±dz = J dx± j ^dz = (J^j J dx± 

where g a b\s * s the induced metric on the trapped surface, which can be calculated by using 
du = 0, dv + dipi(z) = in the line element. 



2. After the collision for weak shocks 

For the line element ( 12. 5 p we parametrize the location of the apparent horizon by r = r^r) 
at constant r, which leads to l^dx^ = (cj — c^r'y) dr + c 2 dr . The conditions for null vectors 
are c\ = —\c 2 A + c 2 r' h and c 2 = 0, so that for constant r 

£ = (-~Ao,o,o,iJ /; = (1,0,0,0,0). 

Vanishing expansion implies E = <9 T £ + |A9 r E = 0, or equivalently S 1 = where we recall 
that 5 1 = E 3 . Because of ( 12.11c) . an equivalent condition is 6{r = r h ) = 0, which is 



sometimes easier to use because the definition of 9 does neither involve the function A nor 
explicit time derivatives. 

Using the values for A, £ from (12. 6p to solve S(r = 77J = 0(fi 3 ) (or equivalently inte- 
grating S to obtain 9{r) and solving for 8(r = r^) = 0(fi 3 )) one finds in the limit r — > 

rfc(r = 0)~ (jjp) ' (2.13) 

and the horizon area becomes 

Mr = 0) = J dx ± d V E 3 ~ ^ V- , (2.14) 

with V = J dx±drj. 

This result is qualitatively the same as in Sec. Ill B 11 which is encouraging. However, 
there may be sizeable quantitative corrections to the above numbers, which can be traced 
back to the approximation used in deriving (12.61) . namely small \i. In terms of the variable 
z = 1/r it is apparent that while corrections C(/i 3 ) to (12. 6 j) are suppressed close to the 
boundary z — 0, they become of order unity when fiz 3 ~ 1, or r ~ /2 1//3 . Hence, the line 
element is not a valid approximation close to the apparent horizon and in particular will 
not fulfill Einstein's equations there. For this reason, we will make an ansatz for the line 
element in Eddington-Finkelstein coordinates that corresponds to (12. 6p for small /i, but is a 
solution to Einstein's equations everywhere. 



C. An ansatz for the post-collision line element 

The ansatz we choose for the line element is to take 

+ (7^^±^, (2 , 5) 
105(r 4 + c/Zs) 

where c is a positive constant ('fudge parameter'). This ansatz agrees with the "perturba- 
tive" result (I2.6P in the limit of small ft and/or small c. In the scheme we will employ, the 
coefficient functions A, B can be calculated numerically from Einstein's equations. We are 
then able to study a toy model for shock wave collisions that involves one unknown number, 
c. 

On the CFT side, the above initial geometry corresponds to a strongly coupled gauge 
theory (N = 4 SYM) with the energy density 

e = T TT = kJi 2 tI + 0(t 3 ) . (2.16) 
8 



TABLE I. Smallest allowed c at different initial times. 



TO 





0.1 


0.2 


0.3 


0.4 


0.5 


1.0 


Cmin(To) 


2.88 


3.43 


4.10 


4.93 


5.97 


7.26 


19.84 




0.83 


0.92 


1.02 


1.15 


1.30 


1.47 


2.85 



This initial condition implies that the initial energy density does not depend on the choice 
for c (nor does any other component of the CFT stress tensor). 

By contrast, the area of the apparent horizon Ah(r) = S(rh,r) depends significantly on 
c. For the gauge choice / = (see (12.91) ) and in the limit tq — > and c< 1, the horizon 
position and area correspond to Eqs. fl 2 . 1 3 [) and (12.141) . respectively, while for c ^> 1 they 
are given by r h {r < 1, c > 1) = (gf^) and A h {r < 1, c > 1) = r 2 h V . Since A h is 
a monotonously decreasing function of c, we may try to approximate its behaviour by the 
Pade inspired ansatz 

A h (r « 1) = fi 2/3 V^= , (2.17) 



648?r 



2X l/3 



fixing the constants k , k\ from the known large and small c limits as k = ^ r ,,,. 

/ 2 N 2/3 

ki = ( ) . We find that this ansatz gives a quite accurate approximation of the 
numerically determined horizon area for arbitrary c in the limit of small To. 

Unfortunately, not all values of c lead to physically acceptable initial conditions, because 
(I2.7ep requires that 3Sdf,S — 2(d r S) 2 < for r > r^, which is not fulfilled for any c. 
Specifically, 3Sd 2 S — 2(d r S) 2 > for r > r v , with ry(r ,c) a monotonously increasing 
function of c. The condition 7v(to,c) < r/ l (r ,c) leads to the requirement c(to) > c m i n (ro), 
with c m i n specified in Tab. [H 



III. NUMERICS 

Near the boundary r — > oo, the metric coefficient functions may be expanded in a power 
series of the following form 

OO OO 7 oo 

A = r 2 V^ B = F-. S=rV^, (3.1) 

n=0 n=0 n=0 

9 



where a = 1, b = — | logr and c = rs, which are determined by the boundary conditions 



2 - 2 - - i 



AU QO = r A , = --logr, = rsr. (3.2) 

Specifically, to the order we will work, we use 

A = A S + A, B = B S + B, £ = £ s + £ , 

9 = 9 s + 9, = s + 0, (3.3) 

where the index s indicates power series expansions. There are only two series coefficients, 



20 



2l|, the 



ai and a 4 , that can not be solved from ( 12 .7p . Via holographic renormalization 
coefficient a 4 is related to the boundary stress tensor by 

T^u = — ^diag ja 4 , a 4 + -<9 r a 4 , a 4 + -<9 T a 4 , -r 2 (a 4 + r<9 T a 4 ) j , (3.4) 

so that we can read off the energy density, longitudinal and transverse pressure of the medium 

as 

3/ta 4 3k . 3« / r \ 

e = — , p L = — (a 4 + rd r a 4 ) , p T = — ^-a 4 - -d T a 4 J . 

In contrast, a x corresponds to the gauge redundancy in (12.91) and, therefore, does not 
appear in any physical quantity. In appendix |A], all the series metric functions A s , B s , S s , 
9 S and (ft s needed in our code are given by taking a\ and a 4 as arbitrary functions of r. The 
series expansions in (13. ip with different gauge choices of a 4 are also related to each other 
according to the transformation in (12. 9 p and (I2.10p with a\ = 2f. 



Numerical method 



Einstein's equations will be solve by the pseudo-spectral method described in Ref. 22] : 
spectral differentiation in r and finite differences in r. The algorithm in the simplest gauge 
choice / = is described in details in Ref. In this case, we would need to impose lower 
cutoff L m j n for the integration domain which needs to fulfill the requirement L m j n < r^(r) for 
all r. We found this approach to work well for late time (near-equilibrium) situations, where 
it is furthermore computationally cheap. However, at early times (far from equilibrium), 
r/j(r) depends strongly on r, and hence it is inconvenient to set up a computational domain 
with fixed L m i n . In these circumstances, one can use an alternative method. Since the 
inside of the horizon is causally disconnected from outside observers on the boundary, the 

10 



computational domain can be chosen to be r > r^r). One can use the diffeomorphism (12. 9p 
to fix rh to a given value, say, unitjo, in which case the only sensible choice for the cutoff 
becomes L min = 1. Besides the lower cutoff it is also necessary to truncate the computational 
domain at large radii at r = L max for numerical reasons discussed in the next subsection. 

As for any pseudo-spectral method, we have to choose the location of grid points (corre- 
sponding to a choice of basis functions), called collocation points 
choose 



ae 



22|. For N + 1 points we 
(3.5) 



where a, b and c are fixed by 



that is, 



ro = L max , r N = L min , and r N / 2 = H, 



H HL max HL m i n -\- L max L mm 

r) TT T T 

£11 J-'max lj min 

b = lo. 



H-L 

tt2 _ : i 

12 Ljnm.a.T. ^ 



mm 
maxLmi n 



(3.6) 
(3.7) 
(3.8) 



Here, H will be chosen to ensure that our algorithm is numerically stable at a relatively large 
time step dr. In the following, we will denote a function / evaluated at any collocation point 
rj by fj. Then the derivative of the function / at r« is given in terms of the differentiation 
matrix D N+1 , 



N 



where 



22 1 



3=0 



D 



N+l,ij 



ab 



■d 



N+l,ij 



(d 



N+l, 



2N 2 +1 



2 ( 1 _ cos 2(^)) 



(-1)' +J 



cj cos(^)-cos( 
2N 2 +1 



i = j = 

< % = j < N 

1 = j = N 



(3.9) 



1 Careful readers will notice that the mass dimension of rh would prohibit us to set it to unity. However, 
one can fix this by introducing an overall dimensionful scale in the problem that will turn out to cancel 
everywhere in physical observables. The definitions given below should be understood in this sense. 
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with Co = c n = 2 and Cj = 1 otherwise. For the sake of numerical accuracy, if 

/( r ) r )lr^oo = /o( r ) r?1 , n > and fo is a function only of r, (3.10) 
we, instead, calculate the derivative of / by 

f: = ^ + r^D N+1 ,^. (3.11) 

Note that since D^+i acting on a constant vector is vanishing, it contains a zero eigenvalue 
and hence is not invertible. For this reason, we consider an alternative version where we drop 
the collocation point j = (corresponding to r = L max , the point closest to the boundary) 
and define anJVxiV matrix 

D N>ij = D N+1>ij , i,j = 1,...N. (3.12) 

Using Dn instead of -Dat+i, we have to supply boundary conditions at r = L max to conserve 
the total number of equations. is invertible and one can numerically solve -D^ 1 , the 
inverse of D^, once for all to save computation time. 

Assuming that A, B and S are known at r, one can first calculate 9 and at r by solving 
(12.1 laft and (12.11b|) with the boundary conditions in (1A5|) and (1A61) . It is of numerical 



advantage to deal with the "residual" metric functions 9 and defined in (I3.3P instead. 
Using the differentiation matrices, the solutions are 

TV 

®3 = ^2 ^NJi^i > 



t=l 

N / N 



3 E Si*8i E D N + i,kB k - <f>' si , (3.13) 



'3 

i=l \ k=0 



where j = 1, 2, • • • , N, 9q = and 0o — 0. The numerical algorithm to solve Einstein's 
equations is then as follows: 

1. Obtain 04, S and B at r + dr by solving the difference equations of (lA4p . (12. 1 lc|) and 
(]2.11dj) . To be more specific, in our code the equations are solved using a third-order 



Adams-Bashforth method, that is, 

dr 

h(r + dr) = h(r) + — [23v(r) - 16u(r - dr) + 5v(r - 2dr)\ , (3.14) 
for a general ordinary differential equation of the form ^ = v. 

12 



2. Then, calculate 64 and ai/r^, at r + dr. To do this, we use the first-order implicit Euler 
scheme to discretize 

d , . , , /i(t + dr) - h(r) d 2 , . . /i(r + dr) - 2/i(r) + /i(r - dr) 

(3.15) 

where h — ax, 64. In this paper, we use the following two algorithms corresponding to 
two different ways to fix the gauge function / = ax/2. 

(a) Alg. I: / = 

In this case, one needs only to solve 64 from the discretized version of the equation 
-Bso| ai= o = an d r h (f = 0) can be calculated by 9(r h ,r + dr) = after one 
has obtained 9 at r + dr in Step 3. 

(b) Alg. II: r h (f) = 1 

In this case, one needs to solve two coupled differential equations given by B s0 = 
B Q and 9n(t + dr) = 9(1, r + dr) = 0. Using the discretization in (13.151) . one 
can express 64 (r + dr) as a function of a±(r + dr) and solve 9n(t + dr) = for 
ax(r + dr) using (ETTBl and flAol) . 

3. Next, calculate 9 and at r + dr by (13.131) with the boundary conditions given by 
9q = 9 s0 and O = (fiso, or equivalently, # = and 4> = 0. 

4. Finally, one can calculate A = A + A s at r + dr by integrating (12.1 lei) . The boundary 
conditions are given by A = A s0 and A' = A' s0 , or equivalently, A = and Aq = 0. 
-Djv+i i s n °t invertible either because it has two eigenvectors with eigenvalue 0. As 
a result, the discretized equation of (I2.11el) gives us only N — 1 linearly independent 
equations, which can be chosen as 

N / N 

E D N+i,iA = mS/Sr 2 + E ^+1,^^ 

i=i V j=o 

-A-A" si , (3.16) 

where % = 2, N and 5" is calculated from (13. lip with n = 3. One needs one more 
equation to solve all A^ with i = 1, 2, ■ ■ • , N, which is given by 

N 

A' = J2 D N+i,oiAi = 0. (3.17) 

i=i 
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A at t + dr can be easily obtained by solving the N linear equations (I3.16p . fl3.17p . 

5. Repeating steps 1-4 one can get the geometry in the bulk at all times. 

To summarize, the numerical algorithm discretizes Einstein's equations using the calcula- 
tional parameters L max , N,dr and fixing either / = 0,L min or r h (f) = 1, L min = 1. The 
continuum Einstein's equations are recovered in the limit L mSuX — > oo,N — > oo,dr — > 0. 

From the description above one can expect that the algorithm in the gauge choice / = 
(Alg. I) should be computationally cheaper than that with r h (f) = 1 (Alg. II). At late times, 
the location of the apparent horizon r^(/ = 0) approaches r = 0. Since L min < r^, this 
implies choosing L min ~ 0. However, this can not always be done. We shall see in the next 
section that for some values of c there are coordinate singularities at r = ry ~ 77, in the 
initial metric set up by f l2.15p . In this case, one has to choose L min > ry and, as a result, to 
stop the code when falls below L m i n . In the following, we will use both algorithms: Alg. 
I for the cases ry <C ^(tq) and Alg. II for the cases ry ~ r^(ro). 



B. Code tests: late time dynamics 

Let us present the performance of the above algorithms in the case where analytic results 
are available: the late time (hydrodynamic) behavior. In this case, the initial conditions for 
the code are given by the following approximate solutions in the gauge a\ — 

r z 3 \ r 



A = r 2 + ^, 5 = --log t + ~ , and S = rr j + r 2 , (3.18) 



4 

where a 4 in the leading-order in - is given by a 4 = —^73 with w a constant. It should be 
emphasized that this initial condition is only an approximate solution to Einstein's equations, 
whereas in the later part of this work we will work with exact solutions as initial conditions. 
Here, we will investigate the time evolution of the error thus made. 

Performing a full-blown numerical stability analysis of our algorithm would be interesting, 
but rather complicated, so we leave it for future work. However, an approximate stability 
criterion can be found by considering the differential operator on the left-hand side of (12.1 left 
and (I2.11dp . that is, d T — ~Ad r . Discretizing this operator we find 

~r \ 5i J + \ dr0 v ) > °a = Ai D N+i,ij- (3.19) 



dr V J 2 
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FIG. 1. The maximum eigenvalue A max in Alg. I for N = 128, L m { n = 0.2 and w r = 4. Left: the 
dependence of A max on H. For L max = 20, 40 and 60, the minimum A max is found at H = 2.7, 2.7 
and 2.6, respectively. Right: dependence on L max with H = 3.0. 

One can argue that the for the algorithm to be stable, the time increment 5t has to be small 
enough that drO^ < Estimating the size of by its maximum eigenvalue A max , we 
find dr < -r-^— . 

As shown in Fig. [TJ the maximum eigenvalue A max depends on the choice of the parameter 
H as well as L max . From this figure, one can see that the choice H = 3 effectively minimizes 
Amax and hence should allow algorithmic stability for larger time increments dr. We adopt 
this choice in following. With A max ~ 5 x 10 3 we therefore expect dr < 10~ 4 to be necessary 
for algorithmic stability. 

In Fig. El we show results by Alg. I with dr — 6.1 x 1CT 5 and Alg. II with dr — 6.1 x 1CT 6 
and 6.1 x 1CT 7 . Using Alg. II with dr = 6.1 x 10 -6 , we get numerically unstable results for 
the horizon position and area. However, by choosing a smaller time step, numerical stable 
results can also be obtained by Alg. II, which agree very well with those by Alg. I. This 
provides evidence for the equivalence of the two algorithms. As a rule of thumb, we find 
that dr — ^ for Alg. I and dr = fo r Alg. II generally ensures numerical stability. 

For numerical reasons, it is difficult to do simulations for L max ^ 100. This can be 
understood from Fig. CEJ where it is shown that X max increases exponentially with L max , 
forcing a similar decrease in the time increment dr. Also, for large L max , one confronts the 
subtraction of two nearly equal numbers in step 3 of the algorithm outlined in Sec. IHI A[ 
Fortunatly, we find in practice that for L maiX w^ 2 T > 10 our numerical results stabilize and 
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FIG. 2. Numerical stability of the two algorithms: gauge choice / = (denoted Alg. I) and gauge 
choice = 1 (Alg. II). Shown are results for the horizon area Ah (left, inset zooms to early 

time behaviour), and the energy density e (right) for L max = 20, N = 128 and H = 3.0. For Alg. I, 
results for one choice of time increment dr are shown (further decreasing dr leaves the result 
unchanged). In order to achieve stable results for Alg. II, we have to decrease dr considerably. 
When this is done, the results from Alg. II match those from Alg. I. 

do not change appreciably when further increasing L max . Thus, we are confident that the 
results reported in the following are close to the continuum limit L max — > oo. Conversely, 
note that for w^ 2 t < 0.1 we would need L max > 100 and therefore cannot report results 
for very early initial times. 

We have also studied the dependence of our results on the number of collocation points N. 
We find that results for N = 64, 128,256 with dr = 1/N 2 are essentially indistinguishable, 
while N = 32 is numerically unstable for dr = 1/N 2 , and differs on the percent level for 
dr = 0.1/N 2 . Thus, we are confident that the choice N = 128 is sufficiently close to the 
continuum N — > oo result and shall adopt this choice in the following. 

Since the initial geometry specified in Eq. (13. 18j) is only an approximate solution to 
Einstein's equations, it is important to check whether time evolution will decrease or increase 
the error. To answer this question quantitatively, we investigate the constraint equation 
(I2.7ep by defining at each r 

5 = max \d 2 r Y> + \B' 2 S I , (3.20) 

M 

where for an exact solution to Einstein's equations 5 = 0. As shown in Fig. [31 5 initially 
is sizeable but decreases as a function of time until eventually stabilizing several orders of 
magnitude below its initial value. This implies that our algorithm approaches the exact 
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FIG. 3. Numerical error 5 in the constraint equation (|2.7e|) as a function of r for N = 64, 128 and 
256. To test the code, we start from initial conditions that do not fulfill Einstein's equations, so 5 
is initially large, but we find that in all cases 5 decreases rapidly as a function of time. 

solution to Einstein's equations as time advances, rather than further deviating from it. A 
physical interpretation is as follows: approximate solutions satisfy Einstein's equations at 
large r, but not close to the horizon position r^. However, in the particular coordinates 
we have chosen, the black hole acts as an absorber of the 'offending' modes, pulling them 
behind the horizon. As a result, the approximate solutions can quickly converge into exact 
solutions. 

The late time hydrodynamic results for e, Ah and are given by 



4 w dro ^ 2^ 3 l + 21og(2) , 

f = — -4- 7/? 

3k r 4 /3 3 r 2 ^ i 8r 8/3 



-3 + 2tt 2 + 24 log(2) - 24 log 2 (2) 

486r 10 / 3 
1 2 + vr + 61o S (2) 

w 



w + 0(t~ 4 ), (3.21) 



rr-l .hydro _ 3 W , 2 + 7T + 6lQg(2) 

h ~ W ~2^ + 24^ 



r 



^-60(-l + log(2) + 121og(2) 2 ) + 187r(l + 61og(2)) , 

2592r 2 1 
w 1 8 + 37r-41og(2) 



hydro 

h -^173-^^ 72^573 

1 /C ttw5 3 25tt 1 7T 2 7 log 2 (2) 

+ „ =7= — + — - - — + — - — — + 



w 2 t 7 / 3 18 3 432 81 7776 162 



TrlogH 21og(w) 251og(2) 
18 27 162 



+ 0(0, (3.23) 
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FIG. 4. Comparison of numerical results from our algorithm ("Spectral"), the algorithm from 
Ref. [8] ("CY") and the 2 nd /3 rd order hydrodynamics. Here, we use tqwI' 2 = 4, L max = 20, 
N = 128, dr = jfa. Shown are r^(/ = 0) (right), energy density e (center) and apparent horizon 
area Ah (left). Performing a least-square fit of the numerical result for e with hydrodynamics we get 
w = 1.0573 u>o- The comparison between numerics and hydrodynamics for Ah and (left, right) 
uses this value of w. 

where r h corresponds to the gauge choice / = and C is Catalan's constant. To test 
the accuracy of our code, we show the comparison between our algorithms, the results from 
Ref. [8] and hydrodynamics in Fig. HI Note that there is very good agreement between the 
algorithm Ref. |8J and the code used in this work. We also find that the energy density 
matches the 3 rd order hydrodynamic result at r > 6. We extract the parameter w governing 
the hydrodynamic behaviour ( I3.2ip by performing a least-square fit to our numerical result 
for e. Using this value of w, Ah matches hydrodynamics at a relatively late time while the 
location of the apparent horizon does so at comparatively earlier times. We recall that the 
initial conditions we had chosen did not fulfill Einstein's equations (see Figj3]), yet at late 
times, we recover hydrodynamics with the correct expansion coefficients (I3.21|3.23|3.23p . 
For exact initial conditions, we expect even better performance. 

To summarize, we have extensively tested our numerical algorithm by studying the de- 
pendence of the results on the numerical parameters L max , N, dr, suggesting that we can 
indeed extract results corresponding to solutions of the continuum Einstein's equations un- 
less starting at very early times r. At late times, our numerical results match the analytically 
known hydrodynamic behaviour as well as those from an independent code. In the following, 
we will now use this algorithm to numerically calculate the solution to Einstein's equations 
2 There is one integral constant £3 in which we can not fix because the 3 rd -order hydrodynamic formula 
of A is still missing in the literatures. 
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for initial conditions modelling the collision of shock waves. 



IV. RESULTS 



In this section, we study the toy model described in Sec. Ill CI using the algorithms in 
the previous section. Using holography, this corresponds to a boost-invariant medium with 
energy density e(r C 1) = k/iVq. The increase of e mimics the 'contracting' stage of 
two nuclei passing through each other in heavy ion collisions. At late times, the system is 
expected to be described by hydrodynamics. Since the system expands along the longitudinal 
direction, one can expect that the medium has to stop contracting, and e should eventually 
decrease in order to match onto the hydrodynamic behaviour. Within our toy model, we 
are able to follow and study all stages of this evolution quantitatively in a strongly coupled 
J\f = 4 SYM medium. 



A. Initial conditions 

Using the ansatz metric function in (12.151) . we set up the initial geometry at time r = r 
by the steps given below. Note that — unlike the test case considered in the previous section 
- the resulting initial condition is an exact solution to Einstein's equations. 



1. Initialize a^. Near the boundary r — > oo, the power series expansion of our ansatz 
metric function £ in (12.151) is the same as that of £ in (12.61) up to O(^j). Therefore, 
CI4 must be given by the same expression as that of the approximate solutions in ( 12. 6ft . 
that is, a 4 (r ) = — 4 ^ 3 T ° . 



2. Then, initialize ai(r ), ai(r — dr) and ai(r — 2dr). Inserting the ansatz in (12. 15[) into 
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(I2.11ap . we can solve 9 analytically in the gauge / = tt = 0, 



e _ r 3 r 4 r //7r(36 + 35^p 5 / 3 T 2 ) 



3 ' 4 10572c 3 / 4 
ft (36^2 + 7c^Y /3 r (18 + Sv^V^)) 
210c 3 / 4 

/2 (36v/2 + 7cVV /3 r (-18 + S^cV^Vs^) 



210c 3 / 4 

v^/x (-36 + 35VB/i 2 / 3 r 2 ) 
420c 3 / 4 



arctan 



arctan 



V2r 



ci/4^1/3 



1 + 



V2r 



log 



cV^V3 

^2/3 + r 2_ y^y/ 3 / 

^2/3 + r 2 + v ^cV^l/3 r 



(4.1) 



Under the transformation (12.91) with / = a\/2 we get #(/) and therefore can obtain 
a\ by solving 9 = 0. In Alg. I, one can skip this step. 

3. Next, initialize 64 (r ), &4(r — dr) and 64 (ro — 2dr). As for a 4 , one can get 64 from the 
power series expansion of B in f 12 . 6 j) . In this way, we have 

2 + 20/x 2 r 6 + rai(4 + ra 1 (3 + ra 1 )) 



12r 4 



4. Initialize -B(t) by integrating fl2.7e[) with boundary conditions given by (12. 6p . 



(4.2) 



5. Finally, obtain 0(to), 0(tq) by (13.131) and A(r ) by solving the linear equations in (13 . 161) 
and (I3~T7D . 

To use the third-order Adams-Bashforth method we also calculate the metric functions at 
ro — dr and To — 2dr by repeating all the steps above. 

As explained in the previous section, it is not possible to choose Ji 1 ^tq = for numerical 
reasons. However, since the initial conditions become less and less reliable for larger r , we 
want to choose Ji 1 ^tq as small as possible such that the numerical algorithm can still be 
applied. The smallest value we achieved in practice was /i 1 / 3 ^ = 0.2. 



B. Transition to hydrodynamic behaviour 

The initial conditions we consider do not exhibit hydrodynamic behaviour at early times. 
This can be clearly seen from the time dependence of the energy density, Eq. fl2.16p . which 
is very different from the hydrodynamic r~ 4 / 3 result. To study the transition from the early 
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FIG. 5. The r— evolution of the energy density and horizon area for Top, 1 / 3 = 0.2 and different 
values of c. Lines labelled 'hydro' represent 3 rd order hydrodynamic results. The w values indicated 
are obtained from hydrodynamic fits to e. 

time behaviour to hydrodynamics, we choose particular values for To, c and then evolve the 
initial conditions in Sec. IIVAI forward in time using our numerical algorithm. 

In Fig. El the apparent horizon area and energy density are shown for initial conditions 
with J^^tq = 0.2 and various values of c. As can be seen from this figure, the energy density 
first increases, reaches a maximum at around p, 1 ' 3 ^ ~ 1, and then starts to decrease. One 
expects the late time dynamics to be described by hydrodynamics, Eq.f l3.2ip . We perform 
a hydrodynamic least-square fit to our results for t to extract the p_ r W at timefl 
p j 1 / 3 r > 7.5. The hydrodynamic results are shown together with the full numerical results 
in Fig. As can be seen, the numerical late-time behaviour of both the energy density as 
well as the horizon area (using the same w values) are very well described by hydrodynamics 
for all chosen values of c. One should note that the energy density is initially independent 
of c, but its late time behaviour differs for different c. This indicates that after some pre- 
equilibrium stage, the system indeed thermalizes, with the overall scale w dependent on the 
non-equilibrium initial conditions. 

The behaviour of the pressure anisotropy is shown in Fig. [HI which for our initial condi- 
tions is Pl/pt = — | at t — 0. One observes that while the system does not exhibit perfect 
isotropy (defined by pr = Pl) for the time extent shown, the pressure anisotropy matches 
the (viscous) hydrodynamic result at around fl^T ~ 3. For all practical purposes, the sys- 
3 We have checked that the extracted values for w change by less than 0.0011% if we perform the fit for 
/2 1 ' 3 t > 6, indicating the insensitivity of the extracted w values. 
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FIG. 6. Left: The r— evolution of the pressure anisotropy, pl/pt for Top, 1 / 3 = 0.2 and various 
values of c. Right: The function F(x)/x for to/U 1 / 3 = 0.2 and various values of c. In both the 
figures, the light grey lines are the 3 rd order hydrodynamic results (hydro). 

tem may therefore be regarded as 'in-equilibrium' for all times thereafter. Conversely, there 
does not seem to be a unique value of Pl/pt above which hydrodynamics is applicable. 
We are also able to directly compare our results with those reported in Ref. . Following 



11). we introduce the quantity = 4jr^ for x = -(— a^) 1 / 4 . which is known within 3 rd 
order hydrodynamics (I3.2ip 



F hydro {x) = 2 1 l-log2 15 -2tt 2 - 45 log 2 + 24 log 2 2 

x 3 9nx 27ir 2 x 2 927n 3 x 3 ' [ ' ' 

As shown in Fig. El our numerical results track the 3 rd -order hydrodynamics solution for 
x > 0.65. This finding is consistent with the result reported in Ref. [ll| for rather different 
initial conditions. 



C. Area scaling and analytic approximations 

While it is numerically hard to send To/i 1//3 — > 0, one can hope to learn about the early 
time behaviour by studying generic values of t$. We thus repeat the above calculations 
for c = 16 for different r , always finding that the late time behaviour is well described by 
hydrodynamics with a parameter w depending on r , c, that is w(t , c). Extracting w(t , c) 
by a hydrodynamic fit to the energy density, we may rescale results for e and using this 
quantity so that the late time behaviour becomes universal. The resulting curves are shown 
in Fig. [3 
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FIG. 7. The dependence of our numerical results on tq with c = 16. Shown are the energy density 
(left) and horizon area (right), scaled by the (fitted) parameter w(tq,c) so that the late time 
behaviour is universal. 

In Fig. [8] we plot the extracted values for w(tq, c) as a function of r . Performing simple 
polynomial fits with degree 2 — 6 we can extrapolate to tq — > 0, finding the value w(0,c = 
16) = 0.885 ±0.02fL 2/9 . This suggests that we may try to obtain an analytical approximation 
to the time dependence of the horizon area Ah as follows: since at tq — > the horizon area 
is given by Eq. f TXTTj) . we find V- l A h {r = 0, c = 16)/w(0,c = 16) 3 ~ 0.61 ± 0.04. This 
certainly seems consistent with Fig. [71 Now knowing the late time behaviour of Ah from 
hydrodynamics and the initial value from Eq. fl2.17l) . we may try to interpolate between these 
two using the ansatz 

T ,_i . / -x / 3 u + Miu;(c)r 2/3 

V A h (r, c) w 3 = v ; 4.4 

1 + diW^CjT 2 ' 6 

where we can fix uo,U\,di by matching the known late and early time behaviour. We find 
Mi = di, di = |(1 — u ), u (c = 16) = 0.61 ± 0.04. The resulting time dependence is close 
to the one found in Fig. [7J although it could be further improved by taking into account 
the known higher order hydrodynamic coefficients. How does the ansatz ( 14. 4p perform for 
different values of c? To this end, let us simply assume that uq = 0.61 ± 0.04 for all values 
of c, that is, the behaviour of the horizon area Eq. (I4.4p would be a universal function. In 
this case, it is easy to predict the value of w from inverting Eq. (I4.4p as 

V^(r = 0,c)^ 1/3 



W(C,T = 0) 



(4.5) 



0.61 ±0.04 

where Ah(ro = 0) is given by Eq. (J27TTJ) . Since we do not have direct access to w(c, r = 0), we 
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FIG. 8. Left: fitted numerical values for w(tq,c) as a function of tq for c = 16 and polynomial 
fit of degree n = 2,3,4,5,6 to the To dependence. Right: extracted numerical values of w(tq,c) 
('numerical') vs. universal prediction (|4.5|) . 

furthermore assume that for all values of c, the ratio w(c, r = 0.2) /w(c, r = 0) = gss^fco 02 ' 
i.e., the same as for c = 16. In Fig. [5] we then compare the predicted universal values for 
w(c, p}' 3 To = 0.2) to the values extracted from our numerical simulations. Surprisingly, the 
predictions from the 'pocket formula' (I4.5p turn out to describe the numerical values almost 
perfectly! Thus it seems that — at least within the class of initial conditions we consider 
- the late time hydrodynamic behaviour is to very good approximation determined by the 
area of the black hole horizon at initial times. 



V. A TOY MODEL FOR THE EARLY TIME EVOLUTION AT RHIC/LHC 

In the preceding sections, we have presented numerical solutions for the time evolution of 
energy density and pressure in a strongly coupled medium that is expanding longitudinally 
in a boost-invariant manner. The initial conditions were chosen such as to mimic those 
following the collision of two shock waves with transverse energy density \x given in Eq. (11. II) . 
where we additionally introduced a 'fudge parameter' c that (together with //) determined 
the area of the trapped surface at r = 0. One may now ask how well these numerical results 
correspond to the experimental situation for heavy-ion collisions encountered at RHIC and 
the LHC Using a = 196, 207 and R = 6.4, 6.6 fm for the atomic number and nuclear radius 
of Au and Pb and a/saT/v = 200, 2760 GeV for the collision energies at RHIC and the LHC 
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we have hrhic — 5.9 GeV 3 , hlhc — 81 GeV 3 . Using our result Eq. (14.51) that relates 

the late time behaviour of the trapped surface to that at early times we can calculate the 

entropy density s = ^A h at 'late' times where hydrodynamics applies as 

2.847Tu 2 / 3 /t 1 / 3 
rVO-97 + c 

In order to interpret this as QCD entropy density at r = 1 fm/c, we first need to fix the 

constant k that is related to the number of degrees of freedom we are simulating. At late 

times Eqns. ( 13. 21113. 231) correspond to s = K7r 4 T 3 with T = ^73- Since it is known that 

this entropy density corresponds to three quarters that of the free case, and knowing that 

s qod = 4( - Nc ~ 1 ^ 7NcNf 7r 2 T 3 we find that we need to set 

_ 4(iV 2 - 1) + 7N c N f 
K ~ 60vr 2 

or k — 0.16 for N c = Nf = 3 in order to model QCD. The corresponding temperatures 
for c — > are then T(r = 1 fm/c) ~ 0.6, 1 GeV for RHIC and LHC energies, which are 
much too large. We may reduce the temperatures to be more in line with values used 



in actual hydrodynamic simulations for RHIC and the LHC (see e.g. Refs. |23H26|) by 
using the fudge parameter c. Additionally, since we are limited to r > 0.2/2 1 / 3 one has to 
rescale the entropy density by a factor of order one (see preceding section). In practice, 
therefore, we choose c = 64 (RHIC) and c = 512 (LHC) which with the measured values 
of wfi~ 2 / 9 = 0.8264, 0.6587 give T(r = 1 fm/c) ~ 0.33, 0.41 GeV for RHIC and LHC, 
respectively. With these parameters, we have a crude, yet fully dynamic model for the 
bulk evolution of the medium created in heavy-ion collisions from r = to the time when 
hydrodynamics becomes applicable. As an example, the evolution of the energy density and 
pressure anisotropy together with the hydrodynamic results are shown in Fig. [9j Note that 
from the pressure anisotropy it seems that hydrodynamics becomes applicable at around 



0.15 fm/c, regardless of the collision energy. A possible application of our result would 
calculation of the non-equilibrium photon/dilepton production along the lines of Refs. 
3o| or Upsilon suppression [27], which we leave for future work. 
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VI. CONCLUSIONS 



In this work, we have provided numerical solutions to a boost-invariant (toy) model of 
shock wave collisions in AdS^, which could be relevant to the problem of heavy-ion collisions 
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FIG. 9. Time evolution of the energy density (left) and pressure anisotropy (right) for RHIC and 
LHC energies, respectively. 

through the AdS/CFT correspondence. Our initial conditions are such that the initial energy 



density evolution is given by the early time analytic solution from Ref . [3| , whereas the early 
time horizon area is controlled by a fudge parameter. Our numerical results indicate that 
the late time energy density behaviour is given by hydrodynamics with a scale parameter 
that is determined by the initial black hole horizon area. More work is needed to decide 
whether this is an artefact of the class of initial conditions we consider or holds true in 
general. Retuning the number of degrees of freedom to make our equation of state QCD- 
like, and freely choosing the fudge parameter we introduced, we are able to provide dynamic 
models for the early time evolution of the bulk medium following heavy-ion collisions at 
RHIC and the LHC, including thermalization of the system. Our results may be useful 
for applications such as calculating non-equilibrium photon/dilepton production and are 
available upon request. 
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Appendix A: Near-boundary behavior of metric coefficient functions 



In this paper, we need to know the following power series expansions of the metric func- 
tions near the boundary r — > oo 

9 1 , 9 . \ a 4 2 + 3r 2 a? + r 3 a? + rai (4 - 9r 4 a 4 ) - 6r 4 (a 4 + 2b 4 ) 
As = r +ai r + -(a 1 -4a 1 ) + - + — 

+ ™\ a (- 32 + r (~37r 2 a 3 - 10r 3 a 4 + 6ai (-12 + 10r 4 (a 4 + 2o 4 ) + 3r 2 a;) + ra\ 
60r t Y 4 

x (-70 + 45r 4 a 4 + 9r 2 a;) + 12r (a[ + r 2 (2a 4 + 4o 4 - 3r6 4 )))) + 243Q 1 7 5 (2512 
+ r (l782r 3 a 4 + 405r 4 af + r 2 a 3 (4289 - 243 (5r 4 a 4 + 3r 2 a;)) - 3ra 2 (-2162 + 15r 2 
x (54r 2 (a 4 + 26 4 ) + 43a; - 3to")) + 2ai (2950 + 9r 2 (-137a; + 3r (r (-36a 4 - 726 4 
+ 5a; 2 + 54r6 4 ) + 5a"))) - 6r (226a; + r (-30a; + r (80a 4 + 1606 4 - 45a' 2 - 138ro 4 
+ 90r 2 6 4 '))))), (Al) 



21ogr 2 1 + rai 4 + 3rai(2 + raQ & 4 1 , , 2 3 

" 3- " 3^ + 18^3 + ^ + I^vF l b4 + T l 5Ur a i 

+ 15r 3 a 4 + 10ra 2 (10 - 3r 2 a;) - 60a! (-2 + 4r 4 6 4 + r 2 a[) + 8 (2r 3 a 4 - 5ra; 

+ 15r 4 6' 4 ))) + 1 6 - (-3712 + r (-1575r 3 a 4 - 405r 4 a^ + r 2 a 3 (-4031 + 1350r 2 a;) 

Zi .LOU/ J 

+ a 2 (-7008r + 45r 3 (l20r 2 6 4 + 73a; - 7ra'[)) - 2a 4 (3742 + 45r 2 (-47a; 
+ r (r (8a 4 + 7a; 2 + 60r6' 4 ) + 7a;'))) + 6r (410a; + r (-70a; + r (-112a 4 
- 1286 4 + 15 (-7a /2 + 2r(6 4 + 7r6 4 / ))))))), (A2) 



and 



1/s 2 + 3rai 1 10 + 9rai -40 - 3rai(20 + 9ra/ 
S s = r 1/3 r + -—77^ tt^- + — — -777^7 + 



+ (A3) 



6r 2/3 9r 5/3 r l62r 8 / 3 r 2 972r n /3 r 3 
where ai, a 4 and 6 4 are functions only of r, which satisfy the following equation 



2 (-2 - 4rai - 3r 2 a 2 - r 3 a 3 + 6r 4 a 4 + 12r 4 6 4 ) 
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From Eqns. fl2.11cj> and (I2.11dp . one can also get 

rr 4 1 1 1 ( \ 1 \ 

S = — + -(2 + 3ra 1 )r 3 + -a 1 (4 + 3ra 1 )r 2 + -a 2 (2 + ra 1 )r + ( — -a?(8 + 3roi) + -ra 4 J 

+ 2 + ra 1 (4 + ra 1 (3 + ra 1 ))-12^ 4 1 ( 3fl 3 _ 18r 4 fl 4 + 15r 2 fl 2 ( _ 1Q + 3r 

30r 4 r 1080r 5 r 2 V 1 1 1 V 

+ ax (-200r + 216r 5 o 4 + 90r 3 ai) - 4 (28 + 3r 2 (-5ai + r 2 (2a 4 + 46 4 + 15r& 4 )))) 

+ 68Q4 Q r 6 r 3 ( 10672 + r (2394r 3 a 4 + 567r 4 a^ + 7r 2 a 3 (1049 - 405r 2 a' 1 ) 
+ 3ra 2 (5078 - 3r 2 (785a; + 21r (36r6 4 - 5a'/))) 

+ 2ai (9490 + 9r 2 (-575a; + 21r (r (4a 4 + 86 4 + 5a; 2 + 30r&' 4 ) + 5a'/))) - 6r (1150a; 
+ r (-210a'/ + t (-304a 4 - 60 86 4 - 3 15a; 2 + 330r6 4 + 630r 2 6 4 '))))) , 

and 



r 3/2 (_ 2 + 3ra 1 ) v ^ ( 12 + ™i(-4 + Zrai))^ 

H ~ o ;n 1" 



3^ 12r 3 / 2 96r 5 /2 

+ 



(-168 - rai(268 + 5rai(38 + 13roi)) + 768r%) (±) 3/2 



384r 7 /2 

+ l —— (5680 + r (2552r 2 a 3 + 771r 3 a 4 + 24ra 2 (259 - 96r 2 a;) 

6144r y / i ' 



5/2 



/IV 

- 32ai (-287 + 144r 2 (2r 2 6 4 + a[)) + 3072 (2r 3 6 4 - ra[ + 3r 4 o 4 ))) ( - J 

+ 1 —tt^ (-226144 + t (-35150r 3 af - 9615r 4 a? + 8r 2 a 3 (-14539 + 7200r 2 a') 

122880T 11 / 2 v v 1 i i V i) 

+ 16ra 2 (-16711 + 480r 2 (l5r 2 6 4 + 17a; - 3ra/)) - 16a! (23083 + 960r 2 
x (-13a; + r (I0r6 4 + 3 (r (af + 5r6' 4 ) + a")))) + 1024r (140a; 

/ lX 7/2 

+ r (-30a" + r (-8a 4 - 546 4 - 45a' 2 + 120r6 4 + 90r 2 o 4 ))))) - . (A6) 
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